H1. Forest structure is negatively correlated with increasing
elevation. H2. Forest structure is negatively correlated with increasing
slope angle. H3. Forest structure is negatively correlated with north
and east aspects. H4. Forest structure is controlled/determined by
topographic position index/ landform classifications.
(forest structural homogeneity is positively correlated with increasing
spatial resolution of the topographic position index) H5. Forest
structure is negatively correlated with increasing Topographic
Roughness. H6. Forest structure will be positively correlated with land
cover and forest type.
The location chosen for this project is Bighorn National Forest in
Wyoming, United States. Bighorn National Forest is approximately 128
meters long and 48 meters wide, encompassing over 404685 hectares of
land (USDA Forest Service 2021). The elevation of Bighorn National
Forest ranges from 1524 meters, to its highest elevation of 4013 meters
at Cloud Peak (USDA Forest Service 2021). Approximately 66% of the land
cover within Bighorn National Forest is forested and made up of
Douglas-fir, ponderosa pine, lodgepole pine, and spruce (Witt 2008). The
sites chosen for this study are located within the National Elevation
Dataset N35 W108, as shown with the black outline in Figure X, which is
approximately 4828 square kilometers.
From August 12-21, 2021, seventy-four plots at 10x10 meters were
selected across the study area of Bighorn National Forest. The sites
chosen were representative of the adjacent area, with their suitability
primarily derived by their elevation to ensure they were within the
montane biogeographic zone (between 1700 and 3000 meters). [Elevations
lower than 1700 meters, and greater than 3000 meters were derived from a
1/3 arc-second digital elevation model (DEM) were contained in polygons
and excluded from site selection consideration. The secondary
consideration for site selection was given to the accessibility of the
site. The majority of the sites chosen were on average 100 meters from a
secondary road (dirt or gravel and not maintained), with few sites being
greater than 100 meters from a primary road (paved and maintained).
Topographic and geographic features were also included in the analysis
of site accessibility, with steep inclines or declines, cliffs, major
water body crossings and unstable or hazardous features generally
excluded from consideration. Sites that had obvious signs of logging
activity or forest fires were also excluded from consideration. The last
criteria for potential sites were that they be analogous to their
surrounding area, and contained at least one mature tree species
(Diameter Breast Height (DBH) must be 10 centimeters or greater). After
the previous criteria were met, final site selection was derived
stochastically by walking a number of paces in one direction determined
by a random number generator. The location produced from the previous
step was used as the Northwest coordinate of the plot.
For each site, a 10m x 10m plot was demarcated. The plots were aligned north-south, east-west, with markers placed at the northwest, northeast, southeast and southwest corners. The preliminary elevation and photos of the plot were taken at the northwest corner (Fig. 3). At each corner, as well as at the center of each plot, a CI-110 Plant Canopy Imager recorded the coordinates, as and a wide-angle image of the plot canopy. Additionally, the Plant Canopy Imager generated canopy data for the plot including the Gap Fraction, Photosynthetically Active Radiation (PAR) Leaf Area Index (LAI), PAR Average, Sunflecks, Canopy Density, and Leaf Angle. Geographic and topographic data about each plot was noted as well as the presence or absence of rocks, boulders, downed woody debris (DWD) or standing dead trees (snags), water and ground junipers. The vegetative species composition of the plot was identified and for trees with a DBH of 10cm or greater, DBH, richness and occurrence was documented.
library(moments)
library(gvlma)
library(raster)
library(lmtest)
library(car)
library(leaps)
library(tidyverse)
library(leaflet)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(broom)
library(sp)
library(sf)
library(vegan)
library(dplyr)
library(spData)
library(stats)
library(evaluate)
library(spatialEco)
library(knitr)
library(insol)
library(tmap)
library(geodiv)
library(GmAMisc)
library(viridis)
library(rgdal)
library(lidR)
library(lattice)
library(rasterVis)
library(RColorBrewer)
library(ggcorrplot)
library(cowplot)
library(googleway)
library(ggrepel)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)
library(citation)
library(MASS)
library(lme4)
knitr::opts_chunk$set(cache=TRUE, warning=FALSE, message=FALSE)
CanopyData <- read.csv("proj_data/CanopyData.csv")
PlotData <- read.csv("proj_data/PlotData.csv")
TreeData <- read.csv("proj_data/TreeData.csv")
TreeHeightsandAge <-read.csv("proj_data/Plot_tree_stats.csv")
DBHsAndDensity <-read.csv("proj_data/DBH_Density.csv")
nf <- read_sf("proj_data/S_USA.AdministrativeForest/S_USA.AdministrativeForest.shp") %>%
st_transform(.,crs=32613)
bighorn <- dplyr::filter(nf,FORESTNAME=="Bighorn National Forest")
DEM <- raster("proj_data/DEM.tif")
Elevation of Bighorn National Forest.
#Read in the LULC Raster
lulc=raster("proj_data/lulc.tif")
# LULC Names
lulc_names =read_csv("proj_data/LULC.csv") %>%
dplyr::select(class=ECOLSYS_LU...2, ID=Value)
# Convert to a Raster
lulc=as.factor(lulc)
# Update the RAT with a Left Join
levels(lulc)=left_join(levels(lulc)[[1]],lulc_names)
#table(values(lulc))
lulc_colors=data.frame(class=levels(lulc)[[1]]$class)
lulc_colors$col=c("#EBA998","#FFFB00", "#4080BF", "#989F60", "#AF5063", "#3200FF", "#A289A8", "#BE20DF", "#971452", "#6250AF" , "#8EB880", "#624074", "#A0B232", "#B232A0", "#FF00E8","#F1197C", "#7C2B4B", "#7C2B72", "#A0D092", "#D095DE", "#70808F", "#10D4EF", "#AF50A7", "#EF103D","#4F30CF", "#FFC300", "#BCFF00", "#00BBFF", "#2080DF", "#DF8E20", "#9A8A86", "#CF3050", "#708F75", "#8F708D", "#95BCDE", "#8F8F70", "#FFCD00", "#EB98E1", "#543E6C", "#DFBA20", "#FF9300", "#00FF43", "#B25932", "#B2327F", "#543B1B", "#1B543B", "#1B3D54", "#805C73", "#B62B85", "#F119B2", "#F18E19", "#591AAF", "#E0FF00", "#E0FF90")
lulcplot<-rasterVis::gplot(lulc)+
geom_raster(aes(fill=as.factor(value)))+
scale_fill_manual(labels=lulc_colors$class,
values=lulc_colors$col,
name="Landcover Type")+
coord_equal()+
theme(legend.position = "bottom", legend.key.size=unit(.5, "line"), legend.text=element_text(size=8))+
guides(fill=guide_legend(ncol=3,byrow=TRUE))
lulcplot
canopy <- CanopyData %>%
group_by(Plot_Number) %>%
dplyr::select(-Sunflecks) %>%
summarise_if(is.numeric, list(mean = function(x) mean(x, na.rm = TRUE),
sd = function(x) sd(x, na.rm = TRUE)))
plots <- PlotData %>%
dplyr::group_by(Plot_Number) %>%
dplyr::summarize(lat=mean(Lat), lon=mean(Lon))
trees<- TreeData %>%
summarize(Plot_Number, Tree_ID, Species_Scientific_Name, DBH, Aprox_Tree_Age) %>%
dplyr:: group_by(Plot_Number)
#GGPlot of Study Site
studysiteplots<-ggplot(plots)+
geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number)) +
scale_color_viridis_c(option="turbo")+
labs(x = "Longitude", y ="Latitude",
title ="Study Site Plots", color ="Plot Number") +
theme(plot.title = element_text(hjust = 0.5))
studysiteplots
#Isolate standard lidar study plots
Lidarplots <-plots[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43,
51,57, 58,59, 65,66, 69),]
#GGPlot of Lidar Plots
lidarplotsgg<-ggplot(Lidarplots)+
geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number))+
scale_color_viridis_c(option="turbo")+
labs(x = "Longitude", y ="Latitude",
title ="Lidar Plots", color ="Plot Number") +
theme(plot.title = element_text(hjust = 0.5))
lidarplotsgg
tree1<-trees %>%
dplyr::select(-Tree_ID, -DBH) %>%
dplyr::group_by(Plot_Number,Species_Scientific_Name) %>%
dplyr::summarise(count=n()) %>%
tidyr::spread(Species_Scientific_Name,count) %>%
dplyr::mutate_all(function(x) ifelse(is.na(x),0,x))
treecounts <-table(trees$Species_Scientific_Name)
view(tree1)
view(treecounts)
tree_coords <- merge(trees, plots)
This graph displays the location of species across plots.
spatial_tc<-ggplot(data=tree_coords, mapping=aes(x=lon, y=lat, color=Species_Scientific_Name))+
geom_point()+
scale_color_viridis_d(option="turbo")+
labs(x = "Latitude", y = "Longitude", title ="Spatial Occurrence of Tree Species", color ="Species Scientific Name") +
theme(plot.title = element_text(hjust = 0.5))
spatial_tc
This bar graph displays the total observations of each species across all plots. The Pinus contorta or Lodgepole Pine was observed the most frequently with more than 350 sightings. The Pinus flexilis, Limber Pine, and the Pinus ponderosa, the Ponderosa Pine each only had one sighting across the plots.
tree_counts_bargraph<-barplot(treecounts,
main= "Species Occurence Throughout Plots",
horiz=TRUE,
xlab="Occurence",
ylab= "Species",
col = c("#DD8317", "#3A9B44", "#47ACC2", "#EACF4F", "#EA5F4f", "#DD1794","#FFC300"),
names.arg=c("Juniperus occidentalis", "Picea engelmannii","Pinus contorta", "Pinus ponderosa", "Pinus flexilis", "Populus tremuloides", "Pseudotsuga menziesii"),
cex.names=0.2, legend=TRUE)
tpi10m <- raster("proj_data/tpi10m.tif")
tpi100m <- raster("proj_data/tpi100m.tif")
tpi1000m <- raster("proj_data/tpi1000m.tif")
TRI <- raster("proj_data/TRI.tif")
slope <- raster("proj_data/slope.tif")
Slope Angle of Bighorn National Forest.
aspect <- raster("proj_data/aspect.tif")
Slope Aspect of Bighorn National Forest.
eastwest <- raster("proj_data/eastwest.tif")
northsouth <- raster("proj_data/northsouth.tif")
SD10m <- sd(tpi10m[],na.rm=T)
SD10m
## [1] 0.4458888
landform_sd10m <- reclassify(tpi10m, matrix(c(-Inf, -SD10m, 1,
-SD10m, -SD10m/2, 2,
-SD10m/2, 0, 3,
0, SD10m/2, 4,
SD10m/2, SD10m, 5,
SD10m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd10m <- as.factor(landform_sd10m)
rat_sd10m <- levels(landform_sd10m)[[1]]
rat_sd10m[["landform_sd10m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd10m) <- rat_sd10m
# Plot the classification
tpi10mplot <-rasterVis::levelplot(landform_sd10m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd10m$landcover,
main = "10M TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd10m[["landform_sd10m"]])))
tpi10mplot
histtpi10m <-hist(tpi10m[])
SD100m <- sd(tpi100m[],na.rm=T)
SD100m
## [1] 2.491715
# Make Landform Classifications
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform_sd100m <- reclassify(tpi100m, matrix(c(-Inf, -SD100m, 1,
-SD100m, -SD100m/2, 2,
-SD100m/2, 0, 3,
0, SD100m/2, 4,
SD100m/2, SD100m, 5,
SD100m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd100m <- as.factor(landform_sd100m)
rat_sd100m <- levels(landform_sd100m)[[1]]
rat_sd100m[["landform_sd100m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd100m) <- rat_sd100m
#writeRaster(landform_sd100m, file="proj_data/landform_sd100m.grd", overwrite=TRUE)
# Plot the classification
tpi100mplot<- rasterVis::levelplot(landform_sd100m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd100m$landcover,
main = "100m TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd100m[["landform_sd100m"]])))
tpi100mplot
histtpi100m <- histtpi100<-hist(tpi100m[])
SD1000m <- sd(tpi1000m[],na.rm=T)
SD1000m
## [1] 29.7036
# Make landform classes
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform_sd1000m <- reclassify(tpi1000m, matrix(c(-Inf, -SD1000m, 1,
-SD1000m, -SD1000m/2, 2,
-SD1000m/2, 0, 3,
0, SD1000m/2, 4,
SD1000m/2, SD1000m, 5,
SD1000m, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform_sd1000m <- as.factor(landform_sd1000m)
rat_sd1000m <- levels(landform_sd1000m)[[1]]
rat_sd1000m[["landform_sd1000m"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform_sd1000m) <- rat_sd1000m
#writeRaster(landform_sd1000m, file="proj_data/landform_sd1000m.grd", overwrite= TRUE) #tif wont work for as.factor
# Plot the classification
tpi1000mplot<- levelplot(landform_sd1000m, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat_sd1000m$landcover,
main = "1000m TPI Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat_sd1000m[["landform_sd1000m"]])))
tpi1000mplot
histtpi1000m<-hist(tpi1000m[])
#### Read in Plots
plotsf<- plots %>%
st_as_sf(coords = c(x="lon", y="lat")) %>%
st_set_crs(.,value=4326)
# Obtain Plot Coordinates in Lidar Projection
lasfile="proj_data/LAS/plot_52_USGS_LPC_WY_Sheridan_2020_D20_13TCK1344.laz"
lproj=readLASheader(lasfile) %>%
crs()
#10m polys
plot_polys <- plotsf %>%
st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs") %>%
st_buffer(dist = 50,endCapStyle="SQUARE") %>%
st_transform(.,crs=st_crs(4326))
# Plot Specific Change
plotgeo <-plot_polys[52,] %>%
st_transform(.,crs=lproj)
# Read in New LAS File
lasPlots <-readLAS(lasfile,
filter = paste("-keep_xy ",paste(st_bbox(plotgeo),
collapse=" "))) %>%
clip_roi(geometry= plotgeo)
plot(lasPlots)